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PREFACE 


This  report  was  prepared  by  Captain  Karen  J.L.  Faran,  Research  and  Development  Coordinator, 
Applied  Research  Branch,  Experimental  Engineering  Division,  U.S.  Army  Cold  Regions  Research 
and  Engineering  Laboratory  .This  work  was  funded  by  the  Directorate  of  Research  and  Development , 
Office  of  the  Chief  of  Engineers,  Project  4A161 1G2  \T24,  Task  FS,  Work  Unit  017,  Winter  Seismic/ 
Acoustic  Energy  Interactions. 

Marie  Burnett,  Institute  of  Geophysics  and  Planetary  Physics,  Scripps  Institute  of  Oceanography, 
La  Jolla,  California,  provided  a  copy  of  the  S  WIS  code  and  some  helpful  notes  on  its  use.  Dr.  Donald 
Albert  and  Randy  McGilvary  of  CRREL  provided  technical  review  of  this  report. 

The  contents  of  this  report  are  not  to  be  used  for  advertising  or  promotional  purposes.  Citation  of 
brand  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial 
products. 
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An  Analysis  of  the  Stress  Wave  in  Solids  (SWIS) 
Finite  Element  Code 

KAREN  J.L.FARAN 


INTRODUCTION 

The  ability  to  analyze  wave  propagation  for  geometrically  complex  circumstances  is  important  in 
calculating  ground  motion  caused  by  earthquakes,  explosions  or  other  sources  of  seismic  waves. 
Analytical  models  derived  using  separation  of  variables  n«:thods  are  limited  in  this  area  because  they 
can  only  solve  problems  with  simple  geometry.  For  more  complex  situations,  it  is  necessary  to  use 
finite  eletTKnt  or  finite  difference  schemes. 

In  1973,  Frazier  (1974)  developed  the  finite  element  code  Stress  Waves  In  Solids,  or  SWIS.  It  has 
been  used  to  solve  several  challenging  problems  because  it  includes  a  variety  of  seismic  propagation 
modes,  including  body  waves,  interface  waves  and  diffraction.  SWIS  is  able  to  simulate  a  number  of 
seismic  pt'.enomena.  Some  examples  are: 

1.  Explosions  in  geologically  complex  formations. 

2.  Spontaneous  earthquake  ruptures  and  near-field  ground  motions. 

3.  Disturbances  in  laterally  varying  earth  models. 

4.  Wave  propagation  through  buried  and  surface  structures. 

SWIS  is  a  versatile  code  in  that  it  can  solve  problems  in  one,  t  xoor  three  spatial  dimensions  in  either 
Cartesian  or  cylindrical  coordinates.  Although  the  code  assumes  lir,;ar  elasticity  and  isotropic 
materials,  it  is  possible  to  solve  problems  in  regions  containing  up  lO  nine  material  types.  The  grid 
generator  has  a  feature  in  which  the  grid  size  may  be  progressively  expanded  at  10%  per  zone  to 
simulate  a  non-reflecting  boundaiy.  Finally,  SWIS  can  solve  f.ither  static,  diffusion  or  wave 
propagation  problems. 

This  report  describes  how  to  use  the  SWIS  code,  which  was  upgraded  at  the  Center  for  Seismic 
Studies  in  1985.  (The  upgrade  was  annotated  in  the  code.)  First,  it  desci'bes  how  to  create  the  input 
file.  A  discussion  of  the  output  files  follow.  Finally,  examples  of  how  SWIS  was  used  to  solve  three 
wave  propagation  problems  are  discussed. 


PROBLEM  INIOALIZATION 

The  numerical  algorithm  in  the  SWIS  code  i  ontains  features  from  both  finite  element  and  finite 
difference  methods.  The  continuum  is  divided,  using  sp:  tial  interpolation  functions  and  a  virtual  work 
pri  nciple,  but  ihe  sequence  is  modeled  after  Langrangia.i  finite  difference  shock  codes.  Also,  theS  WIS 
code  directly  computes  strain  rate,  stress  and  restoring  fv  rces  instead  of  developing  the  conventional 
finite  element  stiffness  matrix. 

To  define  a  stress  wave  problem  for  the  SWIS  code,  the  fcllowing  quantities  are  required  (Frazier 
1974,  ppll-12.- 


1 .  Coordinate  system  designation; 

a.  Number  of  spatial  dimensions  to  appear  in  the  grid. 

b.  Orthogonal  curvilinear  coordinate  system  to  be  employed  in  the  calculations. 

2.  Grid  configuration;  Although  most  grids  can  be  produced  using  the  grid  generator  in  the  code, 
it  is  possible  to  supersede  the  generator  in  local  regions.  Grid  configuration  is  described  by; 

a.  Spatial  location  of  the  node  points. 

b.  Node  map  to  associate  nodes  with  elements 

3.  Boundary  conditions  and  applied  forces;  Each  directional  component  of  each  node  point  is 
assigned  one  of  the  following  constraint  conditions; 

a.  Unconstrained,  with  applied  body  force  or  surface  traction  to  form  an  array  of  nodal  forces. 

b.  Constrained,  with  nodal  displacement  components  constrained  to  follow  a  specified  time 
history. 

4.  Material  properties,  described  by; 

a.  Density. 

b.  Constitutive  properties  (P-wave  and  S-wave  velocities). 

c.  Dimensionless  coefficient  to  regulate  the  damping  of  spurious  high  frequency  numerical 
oscillations. 

5.  Time  stepping  data; 

a.  Start  and  finish  times. 

b.  Time  step,  Ar. 

6.  Starting  condition*:; 

a.  Velocity  and  displacement  with  respect  to  some  reference  frame. 

b.  Stress  at  the  centroid  of  each  element. 

7.  Presentation  of  results; 

a.  Element  and  node  numbers  for  which  results  are  to  be  printed  at  designated  time  intervals. 

b.  Printer  plots  for  displaying  results  at  designated  time  intervals. 

c.  Time  histories  of  individual  node  points. 

d.  Plot  files  producing  graphical  displays  of  the  computed  results. 


nLESUSEDBY  SWIS 

For  both  input  and  output  files,  SWIS  uses  a  two  part  code  for  its  file  names.  The  first  half  is  the 
lettei  “u”  followed  by  a  one  or  two  digit  code  for  the  Fortran  unit  number  used  in  SWIS.  The  second 
half  consists  of  a  two  or  three  letter  description  of  the  contents.  Thus,  file  “u  1  Sin”  is  designated  as  unit 
1 .5  in  SWIS,  and  is  used  as  the  input  file,  and  file  "uShn”  is  the  name  of  unit  8  and  contains  the  time 
history  for  selected  nodes. 

Input  file 

To  run  SWIS,  create  an  ASCII  file,  for  unit  15  titled  ul5in.  The  format  of  this  file  and  variable 
definitions  are  given  in  Appendix  A. 

Output  files 

SWIS  produces  eight  ASCII  files  that  present  computed  displacement  and  velocity  results  in 
different  formats.  By  setting  variables  in  input  file  u  1  Sin  to  appropriate  values,  it  is  possible  to  either 
suppress  printing  or  set  the  time  intervals  for  recording. 

Each  file  can  be  divided  into  several  blocks  of  information.  A  descriptive  summary  and  format 
outline  for  each  of  the  output  files  is  given  below.  Format  A  indicates  a  character  .string,  I  indicates 
an  integer  and  E  represents  exponential  format. 
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1.  u8hn  provides  displacement  time  histones  for  specified  nodes  at  selected  time  intervals. 


Block  I:  Problem  description  (A j. 

Block  2:  Grid  generation  description  (A). 

Block  3:  a.  Number  of  time  stops  (16). 

b.  Number  of  degrees  of  freedom  (16). 

c.  Number  of  nodes  with  recorded  histories  (I6j 

d.  Time  step  (E12.4). 

Block  4:  Node  numbers  for  plot  histoiy  (1 117). 

Block  3:  Displacements  for  each  listed  node,  for  each  time  interval  (8E12.4). 

2.  u9he  contaias  time  histories  of  element  stress  and  displacement.  Block  5  is  printed  for  each 
nth  iteration  (set  in  input  file  ul5in). 

Blo'k  1:  Problem  description  (A). 

Block  2:  Grid  generation  description  (A). 

Block  3:  a.  Number  of  time  steps  (16). 

b.  Number  of  degrees  of  freedom  and  stress  components  (16). 

c.  Number  of  time  histoiy  elements  (16). 

d.  Time  step  (E12.4). 

Block  4:  Element  numbers  for  time  histories  (1 1 17). 

Block  S:  Displacements  and  stress  components  for  e«;h  element  (8Eri;..4). 

3.  ulOg  contains  information  about  the  deformed  grid.  Block  6  is  printed  only  if  a  force  greatc '  than 
O.OCX)!  N  is  applied  to  the  node.  If  time  history  nodes  arc  identified,  both  blocks  7  and  8  are  printed; 
if  no  nodes  are  identified,  only  block  8  is  printed.  Blocks  -12  are  printed  every  nth  iteration  (set  in 
input  file  ulSin). 

Block  1 :  Problem  description  (A). 

Block  2:  Grit*  generation  description  (A). 

Block  3:  a.  Number  of  spatial  dimensions  [ndimtj  (17). 

b.  2'«‘*"'‘(I7). 

c.  Total  number  of  elements  (17). 

d.  Number  cf  diffe  rent  material  types  (17;. 

Block  4:  Coordinates  be  used  in  grid  generation  mapping  (8E12.4). 

Block  5:  a.  P-wave  ve! ,  -  iiv  (El  2.4). 

b.  S-wave  velot  '  (El  2.4). 

c.  Density  (E12.  v. 

Block  6:  a.  Digit  used  to  separate  data  (16). 

b.  Node  coordinates  (8E12.4). 

Block  7:  a.  Digit  used  to  separate  data  (=10)  (16). 

b.  i^ode  coordinates  of  nodes  wit*  *:me  histories  (8E12.4). 

Block  8:  a.  Digit  u.sed  to  separate  data  (=999)  (16). 

b.  Node  coordinates  of  node  1  (8E12.4). 

Bloc  c  9:  Time(E12.4). 

Block  10:  3.2'«”'"‘(I7). 

b.  Material  number  (16). 

Block  1 1 :  Node  coordinate.*;  of  lowest  node  number :»  elements  (8E12.4). 

Block  12:  (Displacer.icnl)+(velocity)*(damping)of lowestiiodemimberinelements(8E12.4). 


4.  ullvn  supplies  data  for  plotting  node  vectors.  Block  5  is  printed  for  every  nth  iteration. 

Block  1:  Problem  description  (A). 

Block  2:  Grid  generation  description  (A). 

Block  3:  a.  Number  of  spatial  dimensions  (17). 
degrees  of  freedom 

c.  Total  number  of  nodes  (17). 

Block  4:  a.  Integer  code  used  for  specifying  nodal  constraints  (Kl. 

b.  Node  coordinates  (3E12.4). 

Block  5;  a.  Time  advance  (E12.4). 

b.  Displacements  and  velocities  (8E12.4). 

5.  ui2  ve  is  supposed  to  provide  data  for  plotting  element  vectors.  Currently,  no  information  is  sent 
to  this  file. 

6.  ul3ln  provides  displacement  and  velocity  information  for  specified  lines  of  nodes.  Block  4  is 
repeated  for  each  line  of  nodes.  Block  5  is  printed  for  eveiy  nth  iteration  of  the  program  (set  in  file 
u  15in).  In  Block  S,  the  items  b,  c  and  d  are  printed  for  eat.')  line  of  nodes.  Furthermore,  displacement 
and  velocities  (item  d)  are  printed  for  each  node  in  the  line. 

Block  i:  Problem  description  (A). 

Block  2;  Grid  generation  description  (A). 

Block  3;  a.  Number  of  dimensions  (17). 

b.  (Number  of  degrees  of  frcedom)*2  (17). 

c.  Number  of  node  lines  (17). 

Block  4:  a.  Node  line  number  (17). 

b.  Number  of  nodes  (17). 

c.  Node  positions  (8E12.4), 

Blocks:  a.  Time  advance  (El  2.4). 

b.  Node  line  number. 

c.  Number  of  nodes  in  line. 

d.  Displacements  and  velocities  for  each  node  (8E12.4). 

7.  uMdiv  provides  the  divergence  and  curl  information  of  the  nodes  specified  in  file  ul31n. 
Information  is  sent  to  ul4div  only  if  information  is  requested  for  lines  of  nodes,  i.e„  if  data  are  sent 
to  file  u  1 31n.  The  output  file  has  only  one  output  fonnat  blo^  k,  which  is  printed  for  each  mh  iteration 
and  for  each  specified  line  of  nodes.  Currently,  ul4div  is  only  printed  for  problems  with  two  spatial 
dimensions  and  with  a  rectangular  mesh. 

Block  1:  a.  Node  line  number  (17). 

b.  Number  of  nodes  in  line  (171. 

c.  Divergence  and  curl  for  each  node  in  the  line  (8E12.41. 

8.  uI6ou!  ..ummarizes  analysis  description,  provides  summ.iry  of  control  parameters,  grid  defini¬ 
tion,  material  definition,  node  constraints  and  output  specifications.  If  so  desired,  u  1 6out  also  contains 
the  computed  results  for  specified  time  intervals.  The  organization  of  this  file  is  self-evident.  An 
example  follows. 
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Example  of  ui6out 

1.  AHALYSIS  DESCRIPTION 
ul51n.ld.2,  ons-dln  prob,  dt-0.01 


2.  CONTROL  PARRHCTERS 


-  Spatial  Representation: 

Number  of  Space  Dimensions  used  .  1 

Number  of  Degrees  of  Freedom  par  Node  .  1 

Number  of  Stress  Components  .  1 

Solution  Coordinate  Designation  .  0 

Order  of  Fourier  Interpolation  in  Axlmuth  ...  0 

-  Time  Control: 

Number  of  Time  Derivlatives  .  2 

Time.  Step  .  0.0100 

Starting  Time  .  0.0000 

Ultimate  Time  .  2.0000 


3.  GRID  DEFINITION 


-  Grid  Generation,  Designator  HAFIZ  ••  2 
Regular  grid,  each  element  O.OS  Mter  long 

grid  sixe:  NEI  100  NEJ  1  NEK  1 
producing:  100  elements  and  101  nodes 

grid  growth  to  element:  IS  0  JS  0  KS  0 

grid  growth  begins  at:  IG  0  JG  0  KG  0 

corner  nodes  of  the  grid  exterior: 

0.00  10.00 


4.  MATERIAL  DEFINITION 


-  Humber  of  Different  Constituents  1 

HAT  DENS  P-VEL  S-VEL  POIS  DAHP 

1  2.7000  6.3000  3.1000  0.3403  0.0000 

-  Haterial  Numbers  Assigned  to  Individual  Elements 
Lines  of  Data  0 


5.  NODE  CONSTRAINTS 


-  Lines  of  Constraint  Data 
NODE  IDNODE  SPECIFIED 
1  1  0.0000 

101  0  1.0000 


2 

CONSTRAINTS 

0.0000  0.0000 

0.0000  0.0000 


0 

0 


0 

0 


6.  OUTPUT  SPECIFICATIONS 


-  Print  Results  at  Interval  .  0 

-  Plot  Deformed  Grid  at  Interval  ....  0 

-  Plot  Node  Vectors  at  Interval  .  0 

-  Plot  Element  Vectors  at  Interval  . .  0 

-  Plot  <  0)  Node  Lines  at  Interval  ..  0 

-  Plot  Time  Histories  of  (  5)  Nodes: 

21  41  61  01  101 


-  Plot  Time  Histories  of  (  0)  Elements; 


NODE 


FTP  BOUND 
0  I  R  TYPE 
R  P  I  1 
C  E  H  12 
E  T  1  2  3 


INITIALIZATION  SUWtARI  R  P  T  H 

-  A  R  Y  A 

SPECIFIED  CONSTRAINTS  NODE  COORDINATES  HASS  NIPT 

SI  S2  S3  Y1  Y2  Y3  E  T 


1100001  0.000 

101  2  0  2  0  0  0  1.000 


0.000  0.000 

0.000  0.000 


0.00 

10.00 


0.00 

0.00 


0.00  0.13  1031 

0.00  0.14 


MOTION  AT  TIME  -  0.0100  (time  step  -  1) 


ODE  BND . DISPLACEMENT . VELOCITY  COMPONENTS 


21 

0 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE4^00 

O.OOOOEfOO 

O.OOOOE+00 

0 

41 

0 

O.OOOOE-fOO 

O.OOOOE+00 

O.000OE«00 

0.0000E4^00 

O.OOOOE+00 

O.OOOOE+00 

0 

61 

0 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE-fOO 

O.OOOOE+00 

0.0000E4^00 

0 

81 

0 

O.OOOOE+00 

O.OOOOE'»00 

o.ooooc«^oo 

0.0000E4^00 

O.OOOOE'fOO 

O.OOOOE+00 

0 

101 

0 

0.7407E-03 

O.OOOOE+OO 

0.0000E400 

0.7407E-01 

O.OOOOEI-OO 

o.ooooe+00 

0 
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SELECTED  EXAMPLES 


To  test  the  SWIS  code,  stress  waves  were  calculated  for  three  wave  propagation  problems:  one¬ 
dimensional  longitudinal  displacement  subjected  to  impulse  loading;  a  cantilever  be':ai  with  an 
impulse  load  iq)plied  along  the  axis,  at  the  unsupported  end;  and  two-dimensional  wave  propagation 
with  a  vertical  impulse  force  (Lamb’s  problem).  The  input  files  and  results  for  these  test  calculations 
follow. 


Example  1:  One-dimeiuaonal  longitudinal  displacement 


Analytical  solution 

The  first  problem  considered  was  that  of  one-dimensional  stress  longitudinal  displacement,  i.e., 
only  displacements  in  thex-direction  were  allowed.  This  situation  ctescribes  the  wave  propagation  in 
the  middle  of  a  large  piece  of  material,  rigidly  constrained  at  one  face  and  with  a  uniform  pressure 
applied  impulsively  at  the  other  (see  Fig.  1 ).  The  material  is  allowed  to  move  only  in  the  direction  of 
the  applied  force,  and  as  a  result,  all  other  displacements  vanish.  The  equations  of  motion,  initial 
conditions  and  boundary  conditions  reduce  to  the  following  one-dimensional  problem: 


_  1  d^u 
dx^  c?  dt^ 


(1) 


initial  conditions: 
boundary  conditions: 


“(«•»)*  I  (■'•»)-» 

tt  (0,  r)  =  0 
/>(f,r)=/>5(t) 


where  u  =  displacement 
t  =time 

X  =  position  along  beam 

c,  =  [(X+|i)/p]  the  longitudinal  or  F-wave  velocity 
P  =  magnitude  of  the  .;onstant  pressure 
5(t)  =  delta  function 
t  =  length  of  the  beam. 


p  =  material  density. 


Figure  I .  Geometry  for  example  1 ,  one-dimensional 
stress,  longitudinal  displacement. 


The  solution  to  1  can  be  found  by  either  a  separation  of  variables  or  by  using  transforms.  The 
latter  technique  gives  the  solutitm  as  (Graff  1975,  pp.  91-94) 
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uiz  = 

H<t 

-H<t- 

'  • '  pc, 

Cl 

Cl 

-H<t- 

(3/-Hjr) 

1 

Cl 

Cl 

1 

f//<r 

~H<t- 

{5/+x:) 

Cl 

Cl 

(2) 


where  //  -a  >  is  the  Heaviside  function,  defined  such  diat 
0,t<a 
I,  t> a 


H <t -a>  - 


Equatimi  2  defines  a  square  wave  [Mopagating  between  the  two  ends  the  mateiial.  widi  wave 

speed  equal  to  the  longitudinal  wave  speed. 


Inpiufile 

The  mesh  created  f(H'thisexanqdewasastringof201  nodes,  lined  in  diex-directioa,  which  crewed 
200  line  elements  (Hg.  2).  Since  displacement  is  restricted  to  only  die  ^'direction,  h  is  uimecessaiy 
to  create  a  three-dimensional  mesh.  If  die  material  is  aluminum,  values  for  dement  lengb,  time  stqi, 
material  properties,  magnitude  of  die  inqiulse  and  dimensions  of  die  f^;ion  are  as  fiiDows: 


time  step  (d/): 
densi^fp): 

F-wave  velocity  (C|): 
5-wave  velocity  (Cj): 
damping: 
inqmlsefmceCfO: 
length  of  re^on(f): 


DT  =  0.005  (ms) 

DENS(l)  =  2.70(Mg/m5) 
VP(l)=6.30(lcm/s) 

VS(I)  =  3.10(km/s) 

DAMP(1)  =  0.0 
VSPEC(2,1)=  1.0(N) 
YGRID(1,2>-YGRID(1.1)  =  10.0(m). 


1  2  3  4  S  6 

§  t  » . 

12  3  4 


201 


200  Pm 


Figure  2.  Finite  element  mesh  for  one-dimensional  stress  prtMem  (200  ele¬ 
ments,  201  nodes). 


For  diis  problem,  node  1  was  assigned  zero  di^lacement  to  meet  the  fixed  end  condition  (line  13 
of  the  following  file).  A  unit  impulsive  force  was  applied  to  the  freeend  of  die  beam,  node  201,  attune 
t = 0  (line  14).  Finally,  recrmls  of  the  di^lacements  were  made  for  five  nodes  along  die  beam:  41, 81, 
121, 161  and  201  (linel7).Theinputfileforthisexamplefollows(entiiesconeqx)odto^)pendixA). 


Entry  Une 
A  1 

B  2 

C  3 

D  4 

E  5 

F  6 


Test  input,  one-dim  prob,  dt=0.005  ms 
11  0  0 

2  0.005  0  10 

Regular  grid,  each  element  0.05  m  long 

200  0  02  2000000 

0.0  10 
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Z-/0 


r=8 


r=2 


b.  SWIS  fdx  =  0.05  m;  dt  =  0.005  s). 
Figure  3.  Comparison  of  analytical  and  SWIS  waveforms  calculated  for  example  1. 


Comparison  of  output  to  theory 

The  disturbance  for  the  given  parameters  should  be  a  square  wave  reflecting  between  the  two  ends 
of  the  material,  at  the  longitudinal  velocity  of  6.3  km/s.  Figure  3  shows  that  the  expected  and  calculated 
waveforms  match. 

One  of  the  shortcomings  of  the  SWIS  solution  is  the  large,  unrealistic  amount  of  ringing  in  the 
results.  Much  of  this  oscillation  has  been  eliminated  hx>m  previous  runs  by  decreasing  both  the  element 
size  and  time  step  (Fig.  4).  It  is  expected  that  the  solution  could  be  further  refined  by  additional 
reductions  in  the  spatial  and  time  increments. 

Another  possible  way  of  reducing  the  oscillations  and  removing  the  high  frequency  noise  in  the 
rigures  would  be  to  introduce  a  damping  factor  with  the  material  parameters.  Figure  5  shows  that  a 
damping  factor  of  0.2  significantly  removes  the  oscillations  in  Figure  3b,  and  the  solution  using  this 
damping  factor  closely  resembles  the  analytical  solution.  This  method  may  have  adverse  effects  on 
the  solution,  however,  in  that  the  higher  damping  factors  change  the  form  of  the  calculated  waves.  As 
seen  in  Figure  3,  the  solutions  obtained  using  non-zero  damping  factors  have  slightly  rounded  comers 
and  finite  rise  times.  However,  the  damping  factors  considered  did  not  seem  to  affect  the  amplitude 
of  the  wave,  nor  did  they  change  the  velocities  at  which  the  disturbances  travel. 
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x=10 


x=8 


1=6 


x=4 


x=2 


4  8 

Time,  ms 


a.  dx  =  0.10  m,  dt  =  0.01  s. 


Figure  4 .  Comparison  of  different  space  and  time  steps for  example  l.The  above  plots 
use  the  same  vertical  scale. 


Example  2:  Cantilever  beam 


Analytical  solution 

The  second  example  considered  was  that  of  a  wave  propagating  along  a  long  and  very  thin  rod,  or 
one-dimensional  stress  where  the  longitudinal  normal  stress  is  a  function  of  position  along  the  rod 
and  time  only  (Fig.  6).  All  other  stresses  vanish,  and  elements  are  allowed  to  deform  in  the  transverse 
direction.  The  equations  of  motion  reduce  to 


d^u  1 


2 

Cb 


dt 


(3) 


9 


Time,  ms 

a.  Damping  =  0.2. 


Time,  ms 

b.  Damping  =  0.4. 


x=10 

x=8 

x=6 

x-4 

1=2 


Figure  5.  Damping  effects  on  one-dimensional  model  (example  1).  The  above  plots 
use  the  same  vertical  scale  (dx  =  0.05;  dt  =  0.005  s). 


Time,  ms 


c.  Damping  =  0.6. 


Figure  6.  Geometry  of  wave  propagationfor  example 
2,  a  cantilevered  beam. 


where  Cb  =  the  beam  velocity 

p 

E  =  |i[(3X+2ii)/(A,+ix)],  Young’s  modulus  of  elasticity 
p  =  material  density. 

With  no  initial  displacements  nor  velocities  along  the  beam,  and  boundary  conditions  of  u(0,r)  = 
0  andP(4r) the  solution  to  this  problem  is  almost  identical  to  the  previous  problem.  The  only 
difference  between  the  two  solutions  is  the  velocity  at  which  the  wave  propagates  through  the  beam 
(C|,  <  c,).  Manipulation  of  the  relations  between  the  material  constants  yield  the  following  relation  for 
C),  in  terms  of  the  longitudinal  and  transverse  velocities 

cb  =  2cf  [(1.5  cf  -  2  cf)  /(c?  -  c,^)] (4) 
Laplace  transform  techniques  (Graff  1975,  pp.  91-94)  give  the  solution  to  the  problem  as 

//</-' - - 1>  - 

Cb 


H<t 


Cb  Cb  f 


H<t 


Cb  Cb 


(5) 


where «  =  displacement 

P  =  magnitude  of  the  load 
p  =  material  density 
Cj,  =  beam  velocity 
i  =  length  of  the  beam 
t  =  time 

X  =  position  along  beam 

H<t-a>  =  Heaviside  step  function,  defined  in  the  previous  example. 

The  solution  to  eq  3  is  a  square  wave  propagating  between  the  ends  of  the  material  at  the  beam 
velocity.  Because  the  beam  velocity  c^,  is  less  than  the  longitudinal  velocity  C] ,  this  wave  travels  slower 
than  the  wave  in  example  1 .  The  amplitude  of  the  resulting  wave,  however,  is  larger  than  that  of  the 
previous  example.  A  plot  of  displacement  versus  time,  for  five  points  on  the  beam,  is  given  in  Figure  8a. 

Input  file 

This  example  differs  from  the  longitudinal  displacement  problem  because  the  nodes  must  be 
allowed  to  move  in  the  transverse  dir^tions  (because  of  the  Poisson  effect).  A  one-dimensional  mesh 
is  not  capable  of  handling  these  displacements,  and  so  either  a  two-  or  three-dimensional  grid  must  be 
used.  To  reduce  computation  time,  a  two-dimensional  mesh  was  created  (Fig.  7)  to  model  an  aluminum 
beam.  The  parameters  for  this  example  follow. 

time  step  (At):  DT  =  0.{X)5  (ms) 

density  (p):  DENS(l)  =  2.70  (Mg/fn^) 

F-wave  velocity  (c , ):  VP(  I )  =  6.30  (km/s) 


II 


5-wave  velocity  (r,): 
beam  velocity  (c,,): 
damping: 
impulse  force  (P): 
length  of  beam  {(): 


VS(l)  =  3.10(km/s) 

C{,  =  5.1  (km/s) 

DAMP(1)  =  0.0 
VSPEC(2,1)=1.0(N) 
YGRID(1,2)-YGRE)(1,1)=  10.0(m). 


{0.+0.05) 
4031 


(0.-0.05) 


(10.+0.05) 


(10,-0.05) 


Figure?. Finite  element  meshofbeamproblem(400  elements, 603  nodes). 
Impulse  force  applied  at  nodes  201, 402  and  603;  nodes  1, 202  and  403 


For  this  example,  the.displacements  for  the  nodes  atx= 0,  nodes  1 , 202  and 403,  were  set  identically 
equal  to  zero  (lines  13, 15  and  17  of  the  following  input  file).  At  the  free  end  of  the  beam,  a  unit  impulse 
was  applied  in  the  direction  of  the  beam  axis  (lines  14, 16  and  18  of  the  following  input  file).  7116 
longitudinal  displacements  were  recorded  for  five  nodes  located  on  the  center  fiber  of  the  beam  (line 
22)  (entries  correspond  to  Appendix  A). 
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A  1  2-d  model  of  beam,  dr=0.005  ms 
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Comparison  of  output  to  theory 

A  plot  of  displacement  versus  time,  as  calculated  by  SWIS,  for  the  above  input  file  is  shown  in 
Figure  8b.  SWIS  calculates  a  waveform  with  a  shape  and  velocity  close  to  that  of  the  analytical 
solution.  As  in  the  case  of  the  one-dimensional  strain  example,  it  is  expected  that  refining  the  input 
mesh  and  decreasing  the  time  step  could  further  improve  the  results. 
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a.  Analytical  solution. 


X=10 
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x=Z.S 
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b.  SWIS  fdx  =  0.05  m;  dt «  0.00S  s). 


Figure  8.  Comparison  of  analytical  and  SWIS  waveforms  calculated for  example  2. 


a.  damping=0.2  b.  dain})ing^0.4 


Figure  9.  Dancing  effects  on  beam.  Plots  use  the  same  vertical 
scale  (dx  *  0.05  m;  dt  =  0.005  s). 


Applying  a  damping  ftetor  again  removes  much  of  thr  oscillations  QFig.  9).  A  fffi:tor  of  0.2, 
however,  already  seems  to  modify  the  solution  in  that  thit  wavefoim  is  no  longer  a  square  wave. 
Increasing  the  damping  factor  from  0.2  to  0.4  removes  more  of  the  high  frequency  components,  but 
results  in  onlyasmallchange  in  thesolution.Forthesedampingfactors,littleornoreductionis  noticed 
in  the  amplitudes  or  the  wave  velocities. 

Example  3:  Lamb’s  problem  in  two<dimensional  Cartesian  coordinates 
The  third  example  treated  the  two^imensional  Lamb’s  problem,  a  vertical  point  load  applied 
impulsively  in  the  plane  of  the  grid  (Fig.  10).  The  results  from  this  example  were  compared  to  the 
waveforms  generated  by  other  computing  schemes. 
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5{t) 


Figure  10.  Example  3,  geometry  of 
Lamb’s  problem. 


Input  file 

No-displacement  constraints  were  applied  to  the  vertical  sides  of  this  mesh  so  that  the  waves  would 
reflect  off  of  the  sides.  To  compare  the  solution  from  SWIS  to  other  models  {Kuhn  1985,  p.  1 1 12),  the 
following  parameters  were  used; 


time  step  (At): 
ending  time: 
density  (p): 

/*-wave  velocity  (cj): 
S-wave  velocity  (c,): 
impulse  force  (F): 


DT  =  2(ms) 

TMAX  =  400(ms) 
DENS  =  1.0  (Mg/m-) 
VP=1.00(km/s) 
VS(l)  =  0.60(km/s) 
VSPEC(2,1)  =  1.0(N). 


For  this  problem,  SWIS  was  run  with  several  input  files  to  observe  the  effect  of  the  damping  factor 
and  to  get  infoimation  for  different  types  of  plots.  For  all  of  the  input  files,  however,  the  mesh  used 
was  a  two-dimensional  grid  with  rectangular  elements,  each  3  by  3  m.  The  mesh  bad  70  elements  in 
each  direction,  and  had  a  total  of  5041  nodes.  A  vertical  force  was  applied  at  the  left  upper  comer  of 
the  mesh  (node  4971)  and  the  vertical  sides  of  the  grid  were  constrained  so  that  these  nodes  had  no 
horizontal  movement  (Fig.  11). 

The  file  shown  below  was  run  to  obtain  information  for  a  contour  plot.  For  every  10  time  steps  (20 
ms),  data  were  recorded  for  15  strings  of  nodes  (line  17  of  the  input  file),  each  string  containing  15 
nodes  (lines  1 8-32).  The  damping  factor  in  this  run  is  0.2  (last  entry  in  line  10)  (entries  conespond  to 
Appendix  A). 


Figure  11.  Finite  element  mesh  used  for  Lamb's  problem. 
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Discussion  of  output 

To  evaluate  the  results,  horizontal  and  vertical  displacements  and  velocities  were  plotted  against 
time  (Fig.  12and  13)andcontourplotsofthedisplacements(Fig.  14)  were  produced.  The  range  and 
depth  scales  for  Figures  12  through  14  were  chosen  to  match  those  of  Kuhn’s  (1985)  figures.  It  is 
important  to  note  that  the  plots  in  Figure  14  may  contain  some  artifacts  attributable  to  the  automatic 
contouring  algorithm.  For  example,  the  contour  plot  of  the  horizontal  displacement  at  80  ms  (Fig. 
14al )  indicates  zero  displacement  at  about  108  m.  This  particular  contour  line  is  not  part  of  the  wave 
front,  but  a  result  of  the  automatic  smoothing  in  the  contouring  algorithm.  Despite  the  artifacts, 
however,  it  is  relatively  easy  to  identify  the  wavefronts  in  the  contour  plots.  The  contour  lines  of 
interest  are  grouped  closely  to  each  other,  and  compose  the  “steep”  portions  of  the  mapping. 

The  plots  of  horizontal  and  vertical  displacement  in  Figure  12  show  the  disturbance  propagating 
through  the  material.  Since  the  wave  velocities  for  the  material  are  known,  it  is  possible  to  determine 
the  arrival  of  each  wave  front.  For  example,  on  the  1  SO-m  trace  of  horizontal  displacements  in  Figure 
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Horizontal 


Vertical 


Figure  13.  Horizontal  and  ycrtical  surface  velocities 
vs  time  for  example  3  (Ax  =  3  m;  dt  =  2  ms;  damping 
=  0.2). 


12a,  a  disturbance  arrives  at  approximately  180  ms.  This  corresponds  to  a  velocity  of  1  km/s,  and 
implies  that  the  disturbance  is  a  pressure  wave.  A  second  wave  front  reaches  the  180-m  range  at 
approximately 300ms,  hasavelocity  of  aboutO.bkrn/s.andcouldbeeithertheshearor  Rayleigh  wave. 
Figure  1 3,  a  plot  of  horizontal  and  vertical  surface  velocities  against  time,  also  shows  the  propagation 
of  the  three  waves.  Finally,  notice  that  the  waves  are  non-dispersive.  This  agrees  with  theory,  since 
the  example  models  a  non-layered  half-space. 

The  Nyquist  frequency,  or  the  highest  ftequency  that  can  be  monitored  owing  to  the  sampling  time 
step,  is  =  l/(2Ar)  =  1/(0.004  s)  =  250  Hz.  From  Figure  12,  the  period  of  the  Rayleigh  wave  is 
approximately  45  ms,andcorrespondstoadominant  frequency  of22Hz.Thisis  an  orderofmagnitude 
smaller  than  the  Nyquist  frequency,  and  so  it  is  reasonable  to  expect  that  the  Rayleigh  wave  is  well 
represented  in  the  plot. 

The  displacement  contours  in  Figure  14  yield  results  consistent  with  theory.  First,  there  are  no 
horizontal  displacements  directly  beneath  the  source  (x  -  0-m  axis),  a  constraint  set  in  the  input  file. 
Disturbances  at  the  depths  of  80  and  1 60  m  are  observed  on  the  = 0  axis  of  the  vertical  displacement 
contour  plots  at  80  and  160  ms  respectively.  These  disturbances  traveled  at  a  rate  of  1  km/s,  and 
probably  correspond  to  the  pressure  wave.  The  second  disturbance,  the  combined  effect  of  the  shear 
and  Rayleigh  waves,  is  observed  near  the  range  of  48  m  on  the  80-ms  plot  and  at  about  96  m  on  the 
1 60-ms  plot.  Finally,  the  displacement  magnitudes,  especially  in  the  horizontal  displacement  contour 
plots,  fall  away  to  zero  with  increase  in  depth  and  indicate  the  presence  of  a  Rayleigh  wave. 

Since  the  compressional  energy  and  shear  energy  are  proportional  to  the  squares  of  the  divergence 
and  curl  of  displacement,  respectively  (Dougherty  and  Stephen  1987,  p.  242),  contour  plots  of  the 
divergence  and  curl  were  created  to  better  observe  the  arrival  of  the  various  wavefronts  at  / = 80  and 
160  ms.  The  equations  used  to  find  divergence  and  curl  are  given  in  Appendix  B. 

The  contour  plots  of  the  divergence  and  curl  facilitate  observation  of  the  wave  fronts.  The 
divergence  contour  plots  show  the  pressure  wave  front  as  being  almost  spherical.  Disturbances  from 
shear  waves  and  surface  waves  are  present  on  the  curl  contour  plots,  but  it  is  difficult  to  distinguish 
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5.  Divergence.  4.  Curl, 

a.  At  time  t  =  80  ms.  Wave  fronts  are  at  the  following  locations:  P-wavc  at  80  m;  S-wave  ai  18  m;  Rayleigh  wave  at  44  m. 

Figure  14.  Contour  plots  of  displacements  and  divergence  and  curl  of  displacements  for  example  3. 
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3.  Divergence.  4.  Curl, 

b.  At  time  t  =  1&)  ms.  Wave  fronts  are  at  the  following  locations:  P-wave  at  160  m:  S-ivave  at  108  m;  Rayleigh  wave  at  88  m. 

Figure  14  (cont'd). 


between  the  two  waves  at  the  sutface  of  the  material  since  the  wave  speeds  are  almost  equal.  At  a  depth 
greater  than  18  m,  however,  the  Rayleigh  wave  displacements  fall  away,  and  only  the  shear  wave 
remains. 

Computation  time 

On  the  ILLIAC  computer,  Frazier  (1974,  p.  65}  estimated  that  the  calculations  were  processed  that 
the  rate  of  0.4  ms  per  two-dimensional  element  per  numerical  time  step.  With 200 time  steps,  and 4900 
elements  (5041  nodes),  each  of  the  mns  took  approximately  30  minutes  of  real  time  on  a  Masscomp 
5550,  a  32-bit  computer  running  at  20  MHz.  At  this  rate,  the  computer  processes  at  approximately  1 .8 
ms  per  element  per  numerical  time  step. 

Dampingfactor 

Frazier,  when  using  SWIS,  used  different  damping  factors  for  the  longitudinal  and  transverse 
waves.  It  is  not  apparent,  however,  how  he  speciried  the  two  factors  in  the  input  file  as  our  version  of 
the  code  does  not  allow  this  option.  At  this  point,  the  magnitude  required  to  reduce  only  the  high 
frequency  noise  resulting  from  numerical  dispersion  has  not  yet  been  determined.  A  value  of  0.2  does 
not  seem  sufficient  because  the  source  wave  oscillates  much  more  than  what  has  been  observed  in  both 
field  work  and  other  mathematical  models.  Damping  factors .  .  to  0.4  and  0.6  reduced  the  amount  of 
oscillation,  but  also  damped  the  results.  Finally ,  a  value  of  0.8  caused  the  disturbance  to  die  out  almost 
immediately. 

Comparison  with  other  models 

Kuhn  (1985)  also  conducted  a  study  of  Lamb’s  problem  in  two  dimensions.  He  used  the  same 
material  parameters  and  numerically  integrated  the  analytical  solution.  In  his  woric,  however,  Kuhn 


Horizontal 


Veiiical 


Figure  15.  Horizontal  and  vertical  ve-  Figure  16.  Vertical  displacement  time  history  calcu- 

locities  (mis)  calculated  by  Kuhn  (after  lated  by  Frazier  (after  Frazier  1974,  p.  67,  his  Fig. 

Kuhn  1985,  p.  1114,  his  Fig.  6a).  4.2),  damping  =  0.8. 
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used  a  different  approximation  for  the 
impulse  force,  solved  the  problem  in  cy¬ 
lindrical  coordinates  and  used  a  mildly 
viscoelastic  material  for  his  half-space, 
lliese  results  are  shown  in  Figure  IS. 

Kuhn’s  figure  shows  surface  velocities, 
and  contains  two  records  of  16  traces 
each,  with  ranges  varying  from  0-180  m. 

The  middle  column  of  numbers  repre¬ 
sents  gain,  which  is  constant  along  each 
trace,  and  allows  the  comparison  of  abso¬ 
lute  amplitudes  between  different  offset 
traces. 

In  all  of  his  calculations,  Kuhn  used 
only  one  source  function,  which  had  a 
dominant  frequency  of  about  20  Hz  (Kuhn  1985,  p.  1 108).  Because  he  solves  Lamb’s  problem  by 
numerically  computing  the  analytic  integral  solution,  his  waveforms  ate  much  cleaner  and  it  is  easier 
to  distinguish  between  the  different  waves.  It  is  difficult  to  see  the  similarities  between  our  solutitm 
(Fig.  13)  and  Kuhn’s  (Fig.  15)  because  the  finite  element  results  contain  much  noise,  resulting  from 
numerical  dispersion.  However,  the  waves  arrive  at  approximately  the  same  time,  and  the  initial  forms 
of  the  waves  are  similar. 

Frazier  ( 1 974,  pp.  65-74)  used  Lamb’s  problem  in  a  two-dimensional  Cartesian  coordinate  system 
to  evaluate  the  SWIS  code  written  for  the  ILLIAC  computer.  As  mentioned  in  the  section  above, 
Frazier  was  able  to  specify  different  damping  factors  for  the  various  waves.  He  also  investigated  the 
effectiveness  of  transmitting  boundary  conditions,  an  option  that  is  not  available  on  our  version  of 
SWIS.  Finally,  Frazier  used  different  parameters  for  his  calculations,  including  a  different  material, 
smaller  time  and  space  steps,  and  a  different  force.  Since  the  parameters  are  so  different  from  those 
in  our  model,  our  comparison  is  limited  to  the  form  of  the  displacements  (Fig.  16). 

As  a  final  comparison,  we  considered  the  calculations  of  Lamb  (Graff  1985,  p.  369).  In  his  analysis 
of  the  half-space  problem,  he  used  a  line  loading  with  a  time  variation  of 

T 

t2+t2 

where  t  is  a  constant.  If  x  is  small,  Z(i)  describes  a  sharp  impulse.  Lamb’s  results  for  the  horizontal 
and  vertical  surface  displacements  from  the  above  loading  ate  shown  in  Figure  17.  The  time  and 
amplitude  scales  are  not  included  in  this  figure,  but  the  first  disturbance  shows  the  arrival  of  aP-wave, 
the  second  corresponds  to  the  5-wave,  and  the  major  response  is  asctibable  to  the  arrival  of  the 
Rayleigh  wave. 


p  S  ruywgh 


Figure  17.  Horizontal  and  vertical  displacements  calcu¬ 
lated  by  Lamb  (after  Graff  1975,  p.  369,  his  Fig.  6.21). 
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APPENDIX  A:  FORMAT  FOR  INPUT  FILE  ulSin 


Listing 

A. HED 

B. NDIMT  NDFNT  MAPXY  NFOUR 

C. NDBYDT  DT  TMIN  TMAX 

D. GRIDH 

£.  NEI  NEJ  NEK  MAPYZ  NBNODES  IS  JS  KS  IG  JG  KG 

F.  ((YGRID(NDIM,NBN)NDIM=1  .NDIMT)  J«N=1  JfflNODES) 

G.  NNCRDC  (if  NNCRDC=0.  go  to  / ) 

^.forNC=ltoNNCRDC: 

NODEC(NC)  (Y(IW,I=1,3)  (DELY(IJ4C),I=1,3)  NANORDCNQ  IANCRD(NC) 
/.  NENNC  (if  NENNC=0,  go  to  K) 

/.forNC=l,NENNC: 

NELN(NC)  (N0DEE(NJ^C)JI=1JWET)  NAEL(NC)  lAEL(NC)  lANE(NC) 

K. NMAT 

L. forN=lJ4MAT: 

MAT  DENS(MAT)  VP(MAT)  VS(MAT)  DAMP(MAT) 

M.  I'^EMATC  (if  NEMATC=0.  go  to  0) 

N. forNC=lJ®MATC: 

NELM(NC)  NEMAT(NC)  NAEMAT(NC)  IAEMAT(NC) 

O.  NNBCC  (if  NNBCC=0.  go  to  Q) 

P.  for  NC=1,  NNBCC  (if  NODEB(NC)>0,  go  to  R) 

NODEB(NC)  NBTYPE(NC)  (VSPEC(I,NC),I=l,3)  NANBC(NC)  IANBC(NC) 

Q. (BCDIR(NCOMP.NAXIS,NC)J^COMPr.l.3)NAXISsl,2) 

R.  INTPRT  INTPG  INTPNV  INTPEV 

S.  NPLTNL  INTPNL  (if  either=0,  go  to  U) 

r.forNL=lJ«>LTNL: 

NDLN(NL)  NANLN(NL)  IANLN(NL) 

U.  NTHPTS,  (NNPRT(I),I=1  .NTHPTS) 

V. NTHELM  NEPRT(I),I=1,NTHELM 

Definition  of  entries 

Entry  formats  are  noted  in  parentheses  (A  =  character  string;  I  =  integer,  and  E  =  exponential 
format). 

A. HED 

(A)  A  character  string  used  to  describe  the  problem;  to  be  used  as  a  heading  on  output.  An  example 
is: 

Input  file  for  uniform  material-lD  with  point  source  at  surface . 

Be  sure  to  leave  a  space  as  the  first  entry  so  the  first  letter  doesn’t  get  readas  acarriage  control  character. 

B. NDIMT  NDFNT  MAPXY  NFOUR 

NDIMT:  (IS)  the  number  of  spatial  dimensions  (1, 2  or  3). 

NDFNT:  (IS)  the  number  of  degrees  of  freedom  per  node. 

MAPXY:  (IS)  designates  the  type  of  spatial  operator;  choices  are  as  follows: 

For  uniform,  rectilinear  grid  in  Cartesian  coordinates,  MAPXY  =  0. 

For  non-uniform,  skewed  grid  in  Cartesian  coordinates,  MAPXY  =  1 . 


For  non-uniform,  skewed  grid  in  Cartesian  coordinates,  and  to  store  stresses  for  non-linear 
constitutive,  MAPXY  =  2. 

For  cylindrical  coordinates  (r,z)  with  harmonic  interpolation  in  azimuth,  MAPXY  =  5. 
NFOUR:  (15)  Fourier  azimuthal  order  in  cylindrical  coordinates. 

If  MAPXY=5,  NFOUR  =  0. 

C. NDBYDT  DT  TMIN  TMAX 

NDBYDT:  (15)  the  number  of  time  derivatives  in  the  partial  differential  equation; 

For  static  problem,  NDB  YDT=0. 

For  diffusion,  NDBYDT=1 . 

For  wave  propagation,  NDBYDT=2. 

DT:  (FI  0.4)  grid  size  in  time.  DT  should  be  less  than  the  space  grid  size  divided  by  the  longitudinal 
wave  (P-wave)  velocity. 

TMIN:  (F10.4)  starting  time. 

TMAX:  (F10.4)  ending  time  (number  of  time  increments  =  TMAX/DT). 

D.  GRIDH 

(A)  A  description  of  the  grid  generation.  As  with  BED,  leave  a  space  for  the  carriage  control  character. 
An  example  for  an  entry  is: 

Regular  grid,  each  element  10  mx  10  m,  7  km  vertical  by  10  km  horizontal. 

E. NEI  NET  NEK  MAPYZ  NBNODES  IS  JS  KS  IG  JG  KG 
NEI:  (15)  number  of  elements  along  the  /-direction  of  a  block  of  elements. 

NEJ:  (15)  number  of  elements  along  the  /-direction  of  a  block  of  elements. 

NEK:  (15)  number  of  elements  along  the  X-direction  of  a  block  of  elements. 

MAPYZ:  (15)  designates  the  mapping  from  the  curvilinear  problem. 

For  identity  mapping,  MAPYZ  =  0.* 

For  bi-quadratic  mapping,  MAPYZ  =  2. 

For  cylindrical  coordinate  mapping,  MAPYZ  =  3. 

For  spherical  coordinate  mapping,  MAPYZ  =  4. 

NBNODES;  (15)  number  of  nodes  that  are  specified  along  exterior  comers  of  the  grid  (NBNODES 
=  4  is  a  typical  entry). 

IS,  JS,  KS:  (15)  starting  numbers  for  expanding  the  grid  size  at  1 0%  per  zone.  Grid  elements  less 
than  IS,  JS  and  KS  are  progressively  expanded. 

IG,  JG,  KG:  (15)  starting  numbers  for  expanding  the  grid  size  at  10%  per  zone.  Grid  elements 
greater  than  IG,  JG  and  KG  are  progressively  expanded. 


F.  ((YGR1D(ND1MMN)NDIM=1  NDIMT)NBN=  I  NBNODES) 

YGRID(NDIM,NBN):  (FI  0.4)  coordinates  of  the  nodes  at  the  exterior  comers  of  the  grid,  specified 
in  the  order: 


ForW*^min-^min-NBN=l. 
^max'  *^min»  ^min»  NBN  =  2. 
^min'  ^min’  NBN  =  3. 

P0'‘W4ax.^mimNBN  =  4. 

For/^i„./„i„,X^„NBN  =  5. 


*  MAPYZ  =  0  is  not  operational;  use  MAPYZ = 2  for  identity  mapping. 
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An*x’  •^min’  ^nuuc*  NBN  =  6. 

^rain'  '4mk*  ^nux’  NBN  =  7. 

^max’  •^nuK*  ^max*  NBN  =  8. 

Supply  NBN’s  depending  on  the  dimensionality  of  the  problem: 

For  a  one-dimensional  problem,  NBN  =  1-2. 

For  a  two-dimensional  problem,  NBN  =  1-4. 

For  a  three-dimensional  problem,  NBN  =  1-8. 

For  a  two-dimensional  problem,  SOOO  units  along  the  top  and  6000  m  deep,  an  entry  could  be  (the 
numbering  used  for  generating  the  grid  need  not  align  with  the  coordinate  axes,  Y1,Y2,Y3): 

0.0  -6000.0  5000.0  -6000.0  0.0  0.0  5000.0  0.0 

For  cylindrical  coordinates,  enter  the  radius  first,  and  then  the  angle  in  radians.  To  specify  a  full  circle 
(2n  radians)  with  radius  of  10,  the  entry  would  be: 

0  0  10  0  0  6.2832  10  6.2832 

G. NNCRDC 

(15)  number  of  lines  (sequences)  of  data  used  to  supersede  none  coordinates.  If  NNCRD = 0,  skip 
to  entry  /. 

H. NODEC(NC)  Y(IJSC)  (DELY(JJ^C).1  =  1,3)  NANCRD(NC)  lANCRD(NC) 

OPTIONAL.  Specify  node  sequence  only  if  NNCRDOO! 

Complete  for  NC  =  1  to  NNCRDC: 

NODEC(NC):  (15)  first  node  number  of  sequence  on  line  NC. 

Y(I,NC),I=  1,3):  (3F10.4)  coordinates  of  node  number  NODEC(NC). 

(DELY(I,NC),I  =  1,3):  (3F10.4)  increment  to  be  added  to  the  node  coordinates  for  generating 
additional  nodes  in  the  sequence. 

NANCREKNC):  (15)  number  of  additional  nodes  in  sequence  NC. 

IANCRD(NC):  (15)  increment  to  be  added  to  the  node  numbers  to  identify  subsequent  nodes  in 
sequence  NC. 

I. NENNC 

(15)  number  of  sequences  (lines)  ofdatausedto  supersede  node  numbers  associated  with  individual 
elements.  SET  NENNC  =  0  and  go  to  entry  K\ 

J. NEmNC)  NODEE(NJ^C)  NAEL(NC)  IAEL(NC)  lANE(NC) 

If  NENNC  =  0,  do  not  enter  values.  Currently,  the  code  only  reads,  and  does  not  process  these 
variables. 

K. NMAT 

(15)  number  of  materials  being  specified.  1  ^  NMAT  ^  9. 

L. MAT  DENS(MAT)  VP(MAT)  VS(MAT)  DAMP(MAT) 

Specify  properties  for  each  material,  N  =  1  to  NMAT. 

MAT:  (15)  material  number,  1  S  MAT  S  9. 

3ENS(MAT):  (F10.4)  mass  density  for  material  number  MAT. 

VP(MAT):  (F10.4)  F-wave  velocity  for  material  number  MAT. 

VS(MAT):  (F10.4)  S-wave  velocity  for  material  number  MAT. 

DAMP(MAT):  (F 1 0.8)  dimensionless  damping  coefficient  to  suppress  high-frequency  amplitudes 
from  numerical  dispersion. 
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M. NEMATC 

(15)  number  of  assignment  sequences;  assigns  material  numbers  to  elements.  NEMATC=0  for  a 
uniform  material.  If  NEMATC  =  0,  skip  to  entry  O. 

N. NELM(NC)  NEMAT(NC)  NAEMAT(NC)  lAEMAT(NC) 

Do  not  enter  values  if  NEMATC  =  0  (uniform  material).  Enter  values  for  NC  =  1  to  NEMATC. 
NELM(NC);  (110)  first  element  in  sequence  NC. 

NEMAT(NC);  (110)  material  number  for  sequence  NC. 

NAEMAT(NC):  (110)  number  of  additional  elements  in  sequence  NC. 

IAEMAT(NC):  (110)  increment  in  element  number  for  identifying  subsequent  elements  in  the 
sequence. 

O. NNBCC 

(15)  number  of  sequences  used  to  constrain  no^s.  The  number  of  different  “boundary  conditions,” 
such  as  applied  forces  or  displacements,  or  both.  If  NNBCC  =  0,  go  to  entry  Q. 

P. NODEB(NC)  NBTYPE(NC)  VSPEC(WC)J=U)  NANBC(NC)  lANBC(NC) 

Used  to  specify  constraints;  enter  values  for  NC=1  to  NNBCC. 

NODEB(NC):  (110)  first  node  in  sequence  NC.  To  apply  a  rotation  to  a  node,  enter  the  negative 
of  the  node  number. 

NBTYPE(NC):  (110)  multi-digit  constraint  code  for  interpreting  components  of  the  values 
specified  by  VSPEC(I,NC).  The  ones  digit  of  NBTYPE  pertains  to  I = NDFNT;  the  tens  digit  pertains 
to  I  =  NDFNT-1 ,  etc.  The  individual  digits  are  interpreted  as  follows: 

0:  VSPEC  is  an  applied  force. 

1:  VSPEC  is  an  applied  displacement. 

Thus,  for  NDFNT  =  2,  NBTYPE  =  00010  indicates: 

VSPEC(  1  ,NC)  =  displacement  assigned  to  component  #1 . 

VSPEC(2,NC)  =  force  applied  to  component  #2. 

Whereas,  for  NDFNT  =  3,  NBTYPE = (XX)10  indicates: 

VSPEC(1,NC)  =  force  applied  to  component  #1. 

VSPEC(2,NC)  =  displacement  assigned  to  component  #2. 

VSPEC(3,NC)  =  force  applied  to  component  #3. 

(VSPEC(I,NC),I  =  1 ,3):  (3F10.4)  the  value  for  the  /th  component  of  the  force  or  displacement  (as 
specified  by  NBTYPE). 

lANBC(NC):  (110)  increment  in  node  number  for  the  subsequent  nodes. 

If  NODEB  >  0,  go  to  entry  R. 

Q.  (BCDIR(NCOMPJ^AXIS^C),  NCOMP=U),  NAXIS^U) 

(F10.4)  used  to  specify  rotations,  but  not  fully  operational;  vectors  to  specify  rotated  directions  for 
degree  of  freedom  NCOMP  with  respect  to  axis  NAXIS. 

R. INTPRT  INTPG  INTPNV  INTPEV 

Used  to  specify  print  control.  Set  the  value  =  0  to  suppress  the  plot. 

INTPRT:  (15)  interval  between  time  steps  for  printing  computed  results  to  unit  16,  file  ‘ul6out’. 
Set  INTPRT  <  0  to  plot  intermediate  values. 

INTPG:,  (15)  interval  between  time  steps  for  plotting  deformed  grid  to  unit  10,  file  ‘ulOg’.  Set 
INTPG  <  0  to  plot  only  the  undeformed  grid. 

INTPNV:  (15)  interval  between  time  steps  for  plotting  node  vectors  to  unit  1 1,  file  ‘ul  Ivn’. 
INTPEV:  (15)  interval  between  time  steps  for  plotting  element  vectors  to  unit  12,  file  ‘ul2ve’. 
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S. NPLTNL  INTPNL 

Plot  along  specified  lines  of  nodes  (sends  output  to  unit  13,  file  ‘ul31n’  and  unit  14,  file  ‘ul4div’). 
NPLTNL:  (15)  number  of  node  lines  (set  NPLTNL  -  0  to  suppress  plots). 

INTPNL;  (15)  interval  between  time  steps  (set  INTPNL  =  0  to  suppress  plots). 

If  either  NPLTNL  or  INTPNL  =  0,  go  to  entiy  U. 

T. NDLN(NL)  NANLN(NL)  IANLN(NL) 

Specify  lines  of  nodes  to  plot  displacement;  enter  values  for  NL  =  1  ,NPLTNL.  Do  not  enter  values 
if  NPLTNL  =  0. 

NDLN(NL);  (110)  first  node  number  in  the  line  NL. 

NANLN(NL):  (110)  number  of  additional  nodes  in  the  line. 
lANLN(NL):  (110)  increment  in  node  number  along  the  line. 

U.  NTHPTS,  (NNPRT(l).  1  =  IMHPTS) 

Used  to  plot  time  histories  of  node  displacement  to  unit  8,  file  ‘u8hn’. 

NTHPTS;  (110)  number  of  nodes  for  which  time  histories  are  to  be  plotted. 

(NNPRT(I),I  =  1, NTHPTS);  (110)  node  number  for  plot  history.  No  entries  are  needed  if 
NTHPTS  =  0. 


V.  NTHELM  (NEPRTd),  1  =  1  J^HELM) 

Used  to  plot  time  histories  of  element  stress  and  displacement  to  unit  9,  file  ‘u9he’. 

NTHELM;  (110)  number  of  elements  for  which  time  histories  are  to  be  plotted. 

(NEPRT(I),I  =  1, NTHELM);  (110)  element  number  for  plot  history.  No  entries  are  needed  if 
NTHELM  =  0. 
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APPENDIX  B:  CALCULATION  OF  DIVERGENCE  AND  CURL 


SWIS  was  modified  so  that  the  divergence  and  curl  would  be  calculated  for  a  two-dimensional 
problem  in  rectangular  coordinates.  The  following  discussion  applies  to  this  specific  case  only. 

For  the  two-dimensional  case,  divergence  and  curl  are  defined  by: 


divergence  (jc)  = 

Bx  By 


curl  (jc)  = 


Bu2  Buj 
Bx  By 


where  x  =  position 

U]  =  displacement  in  jr-direction 
U2  =  displacement  in  y-direction. 


The  divergence  and  curl  were  calculated  using  finite  differences.  The  values  for  the  comer  nodes 
were  calculated  using  forward  differences  for  both  directions;  edge  node  values  resulted  from  a 
forward  diff'jrence  for  the  direction  perpendicular  to  the  edge  and  a  central  difference  along  the  edge; 
and  values  rbr  nodes  in  the  middle  of  the  mesh  were  calculated  using  central  differences  in  both 
directions. 

In  general,  the  forward  and  central  differences  for  a  partial  derivative  are  given  by  (Abramowitz 
andStegun  1972): 


Forward: 


9/0,0  _  /i.o-/o,o 

Bx  h 


+  0(/i2) 


Central: 


9/o,0  _  /l,l  -/-l,  1  +/1.-1  ~/-l.-l 

Bx  Ah 


+  o(;i^) 


where  h  is  the  distance  between  the  sampling  points,  and/ j  is  the  value  of  the  function  at  the  (/th,yth) 
sampling  point.  These  finite  difference  formulas  use  equally  spaced  sampling  points,  as  shown  in 
Figure  B 1 .  For  a  grid  with  non-uniform  spacing,  the  difference  in  coordinates  must  be  used  instead  of 
the  value  h. 


(-1.14 

.(1.1) 

(0.0) 

(1.0) 

(-i.-iT 

*(1.-1) 

a.  Forward  time  dif-  b.Central finite  difference 

ference  sampling  (two  sampling  (four  points), 
points). 


Figure  Bl.  Sampling  points  for  finite  difference  formulae. 
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Nine  different  sets  of  divergence  and  curl  formulae  were  used  forthe  two-dimensional,  rectangular 
mesh.  Each  of  the  four  comer  nodes  required  a  set  of  formulae,  as  did  the  nodes  on  each  of  the  four 
sides  of  the  mesh.  The  final  set  was  written  for  the  nodes  in  the  center  of  the  mesh. 

For  the  following  equations,  variable  definitions  are  given  as: 

nn:  node  number 

net  number  of  elements  in  the  jc-direction  of  the  mesh 

diver(««):  divergence  at  node  m 

curl(nn):  curl  at  node  nn 

disp(i,nn):  displacement  in  the  ith  direction  of  node  nn 

ynode(i,/)«):  coordinate  in  the  ith  direction  of  node  nn 

Node  locations  for  the  following  formulae  are  indicated  on  Figure  B2. 


MI  VIII  IV 


Figure  B2.  Node  locations  for  finite  difference  formulae  (I — 
lower  left  corner;  II — lower  right  corner;  III — upper  left  cor¬ 
ner;  IV— upper  right  corner;  V— lower  edge  of  mesh;  VI— left 
edge  of  mesh;  VII— right  edge  of  mesh;  VIII— upper  edge  of 
mesh;  IX— middle  of  mesh). 


I.  Bottom  left  comer  [««=!]: 

[disp(l,«n+l)-disp(I.nn)]  [di5p(2, nn+nei+\) - disp(2, nn)] 

~  [ynode(  1 ,  nn+\)  -  ynode(  1 ,  nn)]  tynode(2,  nn+nei+ 1 )  -  ynode(2,  nn)] 

curl  {nn)  =  tdisp(2,«n-t-ne/+ 1 ) - disp(2,  nn)]  _  [disp(l,/w+l)-disp(l,ww)] 

[ynode(  1 ,  nn+ 1 )  -  ynode(  1  ,nn)]  [y  node(2,  nn+nei+ 1 )  -  ynode(2,nn)l 

II.  Bottom  right  comer  [nn  =  nei+l]: 

[disp(l ,  nn)  -  disp(l ,  nn-l)]  [disp(2,  wn+nei+1)  -disp(2,  nn)] 
[ynode(l,nn)-ynodc(I,n«-l)]  [ynode(2,nn+«e/+l)-ynode(2,nrt)] 
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curl  („„)  =  tdisp(2,  nn+nei+l)  -  disp(2,  wt)]  [cltsp(l.  nn)  -  disp(l.  wi-l)3 

[ynode(l,  nn)  -  yncKle(l,  /mi-1)]  [ynode(2,  nn+neM)  -  ynode(2,  /ti)] 

III.  Top  left  comer  [nn  =  net  *  (net  + 1)  +  1]: 

diver  (nn)  =  tdisp(l.  nn+l)  -  disp(l>  nn)]  ^  [disp(2.  nn)  - disp(2,  nrt-nei-\)] 
[ynode(l,  /mi+1)  -  ynode(l,  nn)]  [ynode(2,  nn)  -  ynode(2,  nn-nei-1)] 

curl  (nn)  =  [^‘sp(2,  nn)  -  disp(2,  nn-nei-l)]  _  [dispCl,  nn+l)  -  disp(l,  /i/i)] 

[ynode(l,  nn+l)  -  ynode(l.  /mi)]  {ynode{2,  nn)  -  ynode(2,  nn-nei-l)] 

iV.  Top  right  comer  [nn  =  (net  +  l)*nei+  1)]: 

..  ,  .  Idisp(l,/i/i)-disp(l. /MI-1)]  [disp(2,/Mi)-disp(2,/Mi-«c/-l)] 

(flflj  S  I  .  .11.  II  I.  a|»  .  11.—  ^  L.I  I  I— ■■■..!■  . . 

[ynode(l,  nn)  -  ynode(l,  /mi-1)]  [ynode(2,  nn)  -  ynode(2,  nnr-nei-1)] 

curl  (nn)  -  "  <lisp(2,  nn-nei-l)]  _  [disp(l,  nn)  - disp(1 ,  /mi-1)] 

[ynode(l,  nn)  -  ynode(l,  /mi-1)]  Iynode(2,  nn)  -  ynode(2,  nn-nei-l)] 

V.  Nodes  located  on  bottom  edge  of  mesh  [1  <  /mi  <  (nei  +  1)]: 

..  .  .  [disp(l,/Mi+l)-disp(l,/Mi)]  (disp(l,/Mi)-disp(l,/Mi-l)] 

(/J/l  1  S  '  — I— I  ■  ■'■•I  '  I  Ill  ■  I  I 

2*  [ynode(l, nn+l)- ynode(l, nn)]  2*  [ynode(l, nn) - ynode(l, nn~l)] 
[disp(2,  nn+nei+l)  -  disp(2,  nn)] 

^  . .  .  I  . . — . . . 

(ynode(2,  nn+nei+l)  -  ynode(2,  nn)] 

curl  (nn)  -  ~  disp(2, nn)]  [di5p(l,  nn+l)  -  disp(l,  /i/i-l)] 

(ynode(l,  n/i+l) - ynode(l,  n/i-l)]  2*  [ynode(2,  nn+nei+l) -  ynode(2,  nn)] 


curl  (nn)  =  nn+nei+\)- disp(2, nn-nei-\)] _ 2*  [disp( l,nn)- disp(  1 ,  w/i-1 )] 

2*  [ynodcd,  nn) - ynoded.  nn-1)]  [ynode(2.  nn  +nei+l)- ynode(2.  nn-nei-l)] 

VIII.  Nodes  located  on  the  top  edge  of  the  mesh  [nei  *  {nei  +  1)  +  1  <  «n  <  (jiei  +  1)^] 

tdisp(l,  «n+l)  -  disp(l,  nn)\  [disp(l.  nn)-disp(l,  nti-1)] 

2*  [ynode(l,  «m+1)  -  ynode(l,  nn)]  ^  2*  [ynode(l,  nn)  -  ynode(l,  nn-\)] 

^  fdisp(2.  nn)  -  disp(2.  nn-nei-\)] 

[ynode(2,  nn)  -  ynode(2,  nn-nei-\)] 

curl  (nn)  ~ ”H-«e/-l)l  [disp(l, ««+!)- disp(l.  «n-l)] 

cur  nn  [ynode(l,/in+l)-ynode(l,«n-l)]  2*  [ynode(2,  rtn)-ynode(2,  nn-«e/-l)) 

IX.  Nodes  in  the  center  of  the  mesh: 

[disp(l.  nn+nei+2)  -  disp(l.  nn+nei)]  [disp(l.  nn-nei)  -  disp(l ,  nn-nei-2)] 

2*  [ynode(l,  nn+nei+2) -ynod&(l, nn+nei)]  ^  2*  [ynode(l,««-/je/)-ynode(l,  «n-nc/-2)] 

[disp(2,  nn+nei+2)  -  disp(2,  nn-nei)]  [disp(2,  nn+nei)  -  disp(2,  nn-nei-2)] 

*  2*  [ynode(2,  nn+nei+2)  -  ynode(2.  nn-nei)]  2*  [ynode(2,  nn+nei)  -  ynodc(2,  nn-nei-2)] 


curl(nn)  = 


[disp(2,  tm+/ie/+2)  -disp(2,  nn-nei)]  [disp(2,  nn+nei)  -  disp(2,  nn-nei-2)] 

2*  [ynode(l,  nn+nei+2)  -  ynode(I.  nn+nei)]  ^  2*  [ynode(l,  nn-nei)  -  ynode(l,  nn-nei-2)] 


[disp(l,  nn+nei+2)  -  disp(l,  nn+nep]  [disp(l,  nn~nei)  -  disp(l.  nn-nei-2)] 

2*  [ynode(2,  nn+nei+2)  -  ynode(2,  nn-nei)]  ~  2*  [ynode(2,  nn+nei)  -  ynode(2,  nn-nei-2)] 


32 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


reporting  ouroen  tor  this  coiierhion  ot  inrormaucn  is  estimateo  to  avenge  1  hour  per  nsponse,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathenng  a 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestion  for  reduong  this  burden,  to  Washington  Headquarters  Services,  Oirectorale  lor  infomiation  Operations  and  Reports.  t21S  Jefferson  Daws  Highway.  Suite  1204.  Arlington. 
VA  22202-4302.  and  to  the  Olfics  of  Management  and  Budget.  Paperwork  Reduction  Project  <0704-0188).  Washington.  DC  20503. 


1 .  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

November  1991 


3.  REPORT  TYPE  AND  DATES  COVERED 


5.  FUNDING  NUMBERS 


An  Analysis  of  the  Stress  Wave  in  Solids  (SWIS)  Finite  Element  Code 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Cold  Regions  Research  and  Engineering  Laboratoiy 
72  Lyme  Road 

Hanover,  New  Hampshire  03755-1290 


PE:  6.11.02A 
PR:4A!61102AT24 
TA:FS 
WU:017 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

Special  Report  91-21 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

Office  of  the  Chief  of  Engineers 
Washington,  D.C.  20314-1000 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  is  unlimited. 
Available  from  NTIS,  Springfield,  Virginia  22161 


13.  ABSTRACT  {Maximum  200  words) 

The  Stress  Wave  in  Solids  (SWIS)  finite  element  code  is  a  versatile  program  in  that  it  can  solve  problems  in  one,  two  or  three 
spatial  dimensions.  Although  the  code  assumes  linear  elasticity  and  isotropic  materials,  it  can  solve  problems  in  regions 
containing  up  to  nine  different  material  types.  To  demonstrate  its  utility,  SWIS  has  been  used  to  solve  three  classical  wave 
propagation  problems:  one-dimensional  longitudinal  displacement,  impulse  along  the  length  of  a  cantilevered  beam  and 
Lamb’s  problem.  This  report  describes  how  to  use  SWIS  by  summarizing  the  contents  of  the  input  and  output  files. 
Discussions  of  damping  factors,  computation  times  and  comparisons  to  other  solutions  are  also  included. 


14.  SUBJECT  TERMS 

Finite  element  Numerical  modeling  Stress  waves 


IS.  NUMBER  OF  PAGES 


Wave  propagation  he.  price  code 


17,  SECURITY  CUSSIFICATION 
OF  REPORT 

UNCLASSIHED 


18.  SECURITY  Classification 
''*•  I  HIS  PAGE 

L.v  .  .‘.SblFICD 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIHED 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rav-  2-89) 

PtMCflbMbyANSlSM  Z3»-tl 
2M-t02 


